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ABSTRACT 

A  time-accurate  multibody  dynamics  model  of  the  suspension  system  of  a  tracked  vehicle 
is  experimentally  validated  using  a  full-scale  tracked-vehicle  on  an  N-post  motion  simulator.  The 
experiments  consist  of  harmonic  excitations  at  various  amplitudes  and  frequencies  and  ramp 
excitations  of  the  vehicle  road-wheels  (without  the  track),  with  each  road  wheel  under  one  linear 
actuator  of  the  N-post  motion  simulator.  A  high-fidelity  multibody  dynamics  model  of  the  vehicle 
along  with  the  N-post  motion  simulator  is  constructed.  The  multibody  dynamics  model  consists  of 
rigid  bodies,  joints,  rotational  springs  (that  include  non-linear  rotational  stiffness,  damping  and 
friction),  actuators  and  contact  surfaces.  The  rigid  bodies  rotational  equations  of  motion  are 
written  in  a  body-fixed  frame  with  the  total  rigid-body  rotation  matrix  updated  each  time  step 
using  incremental  rotations.  Connection  points  on  the  rigid  bodies  are  used  to  define  joints 
between  the  bodies  including  revolute,  cylindrical,  prismatic  and  bracket  joints.  A  penalty  model 
is  used  to  impose  the  joint  and  normal  contact  constraints.  The  contact  model  detects  contact 
between  discrete  points  on  the  surface  of  each  wheel  (master  contact  surfaces)  and  a  polygonal 
surface  representation  of  the  linear  actuators  dishpans  (slave  contact  surfaces).  A  recursive 
bounding  box/bounding  sphere  contact  search  algorithm  is  used  to  allow  fast  contact  detection. 
An  asperity  friction  model  is  used  for  the  contact  friction  forces.  The  governing  equations  of 
motion  are  solved  along  with  joint/constraint  equations  using  a  time-accurate  explicit  solution 
procedure.  The  time-histories  of  the  suspension  system  rotational  deflections  are  experimentally 
measured  for  the  various  input  motion  excitations.  The  experimental  measurements  are  compared 
with  the  results  predicted  using  the  computational  model.  The  comparison  shows  that  the  model 
can  predict  with  reasonably  good  accuracy  the  test  tracked-vehicle ’s  dynamic  response. 


1.  INTRODUCTION 

Tracked  vehicles  are  used  in  civilian  and  military  off-road 
ground  vehicles.  Tracks  have  better  mobility  characteristics 
than  tires  on  off-road  and  slippery  terrains.  This  includes: 

•  Better  traction  over  firm  terrains  (such  as  pavement) 
that  is  highly  sloped  or  slippery  (e.g.  wet/icy  pavement). 

•  Better  traction  and  mobility  over  soft  terrains  such  as 
snow,  sand,  and  mud. 

•  Better  positive/negative  obstacle  crossing  capability. 
Tracked  vehicles  can  go  over  larger  obstacles  than 
equivalent  wheeled  vehicles.  Typical  obstacles  that 


ground  vehicles  are  tested  against  include  semi-circular 
and  trapezoidal  bumps/potholes. 

Tracks  can  be  used  on  vehicles  of  various  sizes  ranging  from 
large  tractors  and  excavator  to  small  unmanned  ground 
vehicles.  Tracks  can  also  be  used  in  conjunction  with  tires  in 
order  to  improve  the  vehicle’s  mobility.  Tracks  can  be 
divided  into  two  types: 

•  Continuous  belt  track.  The  cross-section  of  a  continuous 
belt  track  is  similar  in  construction  to  a  tire.  It  consists 
of  a  rubber  matrix  reinforced  with  steel,  Kevlar  and/or 
polyester  wire/ply  along  the  length  and  width  of  the 
track. 
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•  Segmented  track.  The  track  consists  of  rigid  or  flexible 
segments  connected  using  re  volute  joints.  For  rigid 
segmented  tracks,  the  segments  are  made  of  steel  with  a 
rubber  layer  on  the  track-road  contact  surfaces.  Flexible 
segmented  tracks  have  both  reinforced  rubber  and  steel 
segments.  The  rubber  segments  are  used  for  contact 
between  the  track  and  the  terrain.  The  steel  segments  are 
used  to  connect  the  rubber  segments. 

Tracked  vehicles  are  flexible  multibody  systems  involving 
frictional  contact,  large  rotations  and  large  deflections  (of 
the  flexible  track).  A  review  of  flexible  multibody  dynamics 
modeling  techniques  including  deformation  reference 
frames,  treatment  of  large  rotations,  discretization 
techniques,  finite  elements,  constraint  and  contact  modeling, 
and  solution  techniques  is  presented  in  [1].  A  review  of 
tracked  vehicles  models  was  presented  in  [2].  A  finite 
element  multibody  dynamics  model  for  tracked  vehicles, 
that  can  handle  both  flexible  belt-type  tracks  and  rigid 
segmented  tracks,  was  presented  in  [2].  The  model  uses  the 
following  computational  techniques/elements: 

•  Explicit  time-integration  solver  accurate  for  long 
simulation  times  [3]. 

•  The  algorithm  that  is  used  to  account  for  a  rigid  body 
rotational  motion  [4]  uses  the  total  rotation  matrix 
relative  to  the  inertial  frame  to  measure  the  rotation  the 
rigid  bodies.  This  avoids  singularity  problems 
associated  with  3  and  4  parameter  rotation  measures. 
Rigid  body  rotational  equations  of  motion  are  written  in 
a  body  (material)  frame.  Thus,  the  inertia  tensor  is 
constant.  The  equations  of  motion  are  integrated  to  yield 
the  incremental  rotations  angles.  The  total  body  rotation 
matrix  is  updated  using  the  rotation  matrix 
corresponding  to  the  incremental  rotations  angles. 

•  Penalty  formulation  for  modeling  joint  constraints 
including  spherical,  revolute,  cylindrical  and  prismatic 
joints  [4]. 

•  Normal  contact  modeled  using  a  penalty  formulation  [5, 

6]. 

•  Frictional  contact  modeled  using  an  accurate  and 
efficient  asperity-based  friction  model  [7]. 

•  General  contact  search  algorithm  for  finding  the  contact 
penetration  between  finite  elements  and  other  elements 
as  well  as  general  triangle  and  quadrilateral  surfaces  [2]. 

•  Total-Lagrangian  spring,  truss,  beam  and  solid 
nonlinear  finite  elements  with  Cartesian  coordinate 
degrees  of  freedom  [3,  8-10].  This  allows  arbitrarily 
large  element  rotations  with  no  solution  drift  due  to  the 


use  of  displacement  increments.  The  element  DOFs  are 
referenced  to  a  global  inertial  reference  frame. 


•  The  thin  beam  elements  use  a  torsional- spring 
formulation  with  3 -nodes  per  element  [3,  8]  for 
modeling  the  bending  behavior  of  the  beams. 

•  8 -node  natural  deformation  modes  brick  solid  elements 
are  used  [9,  10].  Those  elements  can  also  be  used  to 
model  shells  and  beams.  One  element  through  the 
thickness  is  sufficient  to  accurately  model  the 
membrane,  shear,  and  bending  characteristics.  The  brick 
elements  do  not  exhibit  locking  or  spurious  modes 
(widely  used  techniques  to  alleviate  locking  such  as 
hourglass  control  lead  to  elements  that  do  not  maintain 
solution  accuracy  over  very  long  solution  times).  Any 
material  law  can  be  used  with  those  elements  including: 
linear  elastic,  hyper-elastic,  and  non-linear  laws.  Those 
elements  are  used  to  model  the  rubber  matrix  of  the 
continuous  belt-type  tracks.  The  track  reinforcements 
for  those  tracks  can  be  modeled  using  the  thin  beam 
torsional-spring  element. 


•  Finear/rotational  suspension  springs  that  can  include  a 
prescribed  non-linear  force-deflection  relation,  a  force- 
deflection  rate  relation  and  a  friction  coefficient- 
deflection  rate  relation. 


Figure  1.  Tracked  vehicle  on  the  n-post  motion  simulator. 


The  purpose  of  this  paper  is  to  partially  experimentally 
validate  the  tracked  vehicle  model  presented  in  [2].  The 
experimental  validation  was  carried  out  using  a  full-scale 
tracked  vehicle  that  was  instrumented  and  mounted,  with  the 
track  removed,  on  an  n-post  motion  simulator  (Figure  1). 
The  n-post  motion  simulator  consists  of  hydraulic  linear 
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actuators  under  each  road  wheel  that  can  be  independently 
controlled  to  simulate  typical  road  maneuvers  that  the 
vehicle  is  subjected  to  such  as  going  over  bumps/potholes. 
The  time-histories  of  the  angular  deflections  of  the  road 
arms  are  measured  for  various  input  motion  excitations.  The 
experiments  were  used  to  tune  the  force  -  angular  velocity 
relation  (damping)  of  the  road  arms  torsion  bars.  The  results 
of  the  experiments  are  then  compared  with  the  results 
predicted  using  the  computational  model. 

In  addition,  in  this  paper  typical  simulations  of  tracked 
vehicles  that  can  be  performed  using  the  present  model  are 
also  presented.  Also,  a  lumped  parameters  thick  spatial  beam 
element  is  presented.  The  element  has  two  rigid  body  type 
nodes  with  each  node  having  3  translational  DOFs  and  the 
rotation  matrix  is  used  to  measure  the  spatial  rotation  of  the 
node.  The  element  models  each  response  component 
independently  using  simpler  sub-elements  that  include:  three 
mutually  orthogonal  linear  spring  sub-elements  for  modeling 
one  axial  and  two  shear  deformations;  three  mutually 
orthogonal  torsional  spring  elements  for  modeling  one 
torsional  and  two  bending  deformations  of  the  element.  This 
beam  element  can  be  used  to  model  both  segmented  and 
continuous  tracks. 

The  rest  of  this  paper  is  organized  as  follows.  In  Section  2 
the  translational  and  rotational  semi-discrete  equations  of 
motion  are  presented.  In  Section  3  the  finite  element 
formulation  for  the  truss,  thin  beam,  thick  beam  and  solid 
brick  elements  are  presented.  In  Section  4  the  frictional 
contact  techniques,  including  penalty  formulation  for  normal 
contact,  asperity  friction  model  and  contact  search 
algorithm,  are  presented.  In  Section  5  the  penalty  algorithm 
for  imposing  joint  constraints  is  presented.  In  Section  6,  the 
overall  explicit  solution  algorithm  is  outlined.  In  Section  7, 
the  experimental  setup  is  presented.  In  Section  8  the 
comparison  between  the  experiment  and  model  is  presented. 
In  Section  9  typical  numerical  simulations  of  tracked 
vehicles  are  presented.  Finally,  in  Section  10  concluding 
remarks  are  offered. 

2.  EQUATIONS  OF  MOTION 

In  the  subsequent  equations  the  following  conventions  will 
be  used: 

•  The  indicial  notation  is  used. 

•  The  Einstein  summation  convention  is  used  for  repeated 
subscript  indices  unless  otherwise  noted. 

•  Upper  case  subscript  indices  denote  node  numbers. 

•  Lower  case  subscript  indices  denote  vector  component 
number. 

•  The  superscript  denotes  time. 

•  A  superposed  dot  denotes  a  time  derivative. 


Two  types  of  finite  element  nodes  are  used:  point  particle 
nodes  and  rigid  body  nodes.  Point  particle  nodes  have  3 
translational  DOFs  (Degrees  of  freedom).  The  algorithm  for 
writing  and  integrating  the  equations  of  motion  for  spatial 
rigid  bodies  using  an  explicit  finite  element  code  was 
presented  in  [4].  In  this  algorithm,  a  rigid  body  is  modeled 
using  a  finite  element  node  located  at  its  center  of  mass.  The 
node  has  3  translational  DOFs  defined  with  respect  to  the 
global  inertial  reference  frame  and  a  rotation  matrix  defined 
also  with  respect  to  the  global  inertial  frame.  The  use  of  the 
total  body  rotation  matrix  to  measure  rigid  body  rotations 
avoids  singularity  problems  associated  with  3  and  4 
parameter  rotation  measures  [3]. 

The  translational  equations  of  motion  for  the  nodes  are 
written  with  respect  to  the  global  inertial  reference  frame 
and  are  obtained  by  assembling  the  individual  node 
equations.  The  equations  can  be  written  as: 

MKx‘Ki=Fs‘B+Fa‘Ki  (1) 

where  t  is  the  running  time,  K  is  the  global  node  number  (no 
summation  over  K\  K=\->N  where  N  is  the  total  number  of 
nodes),  i  is  the  coordinate  number  (z=l,2,3),  a  superposed 
dot  indicates  a  time  derivative,  MK  is  the  lumped  mass  of 
node  K ,  x  is  the  vector  of  nodal  Cartesian  coordinates  with 
respect  to  the  global  inertial  reference  frame,  and  x  is  the 
vector  of  nodal  accelerations  with  respect  to  the  global 
inertial  reference  frame,  Fs  is  the  vector  of  internal  structural 
forces,  and  Fa  is  the  vector  of  externally  applied  forces, 
which  include  surface  forces  and  body  forces. 

For  each  node  representing  a  rigid  body,  a  body-fixed 
material  frame  is  defined.  The  origin  of  the  body  frame  is 
located  at  the  node  that  is  also  the  body’s  center  of  mass. 
The  mass  of  the  body  is  concentrated  at  that  node  and  the 
inertia  of  the  body  is  given  by  the  inertia  tensor  defined  with 
respect  to  the  body  frame.  The  orientation  of  the  body- frame 
is  given  by  R £  which  is  the  rotation  matrix  relative  to  the 
global  inertial  frame  at  time  t0.  The  rotational  equations  of 
motions  are  written  for  each  node  with  respect  to  its’  body- 
fixed  material  frames  as: 

4,4  =  T 4  +  K,  -  (4  x  (4,4  )L  (2) 

where  IK  is  the  inertia  tensor  of  rigid  body  K ,  Q  and  Q  are 

the  angular  acceleration  and  velocity  vectors’  components 
for  rigid  body  K  relative  to  its  material  frame  in  direction  j 
(j=  1,2,3),  TsKi  are  the  components  of  the  vector  of  internal 
torque  at  node  K  in  direction  /,  and  TaKi  are  the  components 
of  the  vector  of  applied  torque.  The  summation  convention 
is  used  only  for  the  lower  case  indices  i  and  j.  Since,  the 
rigid  body  rotational  equations  of  motion  are  written  in  a 
body  (material)  frame,  thus,  the  inertia  tensor  IK  is  constant. 
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The  trapezoidal  rule  is  used  as  the  time  integration  formula 
for  solving  equation  (1)  for  the  global  nodal  positions  x: 

xkj  =  x KjAt  +  0.5  (4  +  xfKjAt )  (3  a) 

xkj  =  xKjAt  +  0-5  At  +  xfKjAt )  (3b) 

where  At  is  the  time  step.  The  trapezoidal  rule  is  also  used  as 
the  time  integration  formula  for  the  nodal  rotation 
increments: 


d'Kj=d‘-;‘+0.5At(d'Kj+9‘-;‘)  (4a) 

A0lKj=O.5At(0lKj+0‘:Al)  (4b) 

where  A6Kj  are  the  incremental  rotation  angles  around  the 
three  body  axes  for  body  K.  Thus,  the  rotational  equations  of 
motion  are  integrated  to  yield  the  incremental  rotation 
angles.  The  rotation  matrix  of  body  K  ( RK )  is  updated  using 
the  rotation  matrix  corresponding  to  the  incremental  rotation 
angles: 

R'K  =  R'KAI  R(A0'la)  (5) 


where  R{AOlKi)  is  the  rotation  matrix  corresponding  to  the 
incremental  rotation  angles  from  equation  (4b). 

The  explicit  solution  procedure  used  for  solving  equations 
(1-5)  along  with  constraint  equations  is  presented  in  Section 
7.  The  constraint  equations  are  generally  algebraic 
equations,  which  describe  the  position  or  velocity  of  some  of 
the  nodes.  They  include: 


•  Contact/impact  constraints  (Section  5): 


o 

(6) 

Joint  constraints  (Section  6): 

/(M)  =  0 

(7) 

Prescribed  motion  constraints: 

f({x},t)  =  o 

(8) 

3.  Finite  Elements 

3.1  Truss/Spring  Element 

The  truss  element  connects  two  nodes.  The  internal  force 
in  a  truss  element  is  given  by: 


F 


EAn  ix  CAi 

t{\  t{\ 


(9) 


where  E  is  the  Young’s  modulus,  C  is  the  damping  modulus, 
A  is  the  cross-sectional  effective  area,  /  is  the  current  length 
of  the  truss,  /0  is  the  un-stretched  length  of  the  truss. 


3.2  Torsional-Spring  Thin  Spatial  Beam  Element 

The  torsional-spring  beam  element  developed  in  [8]  is 
used  for  modeling  thin  beams.  The  element  has  3  point  mass 
type  nodes  (nodes  which  have  only  translational  DOFs).  A 
beam  element  is  shown  in  figure  2a  The  beam  element 
connects  the  point  px  (mid-point  of  12)  to  point  p2  (mid-point 
of  23).  The  slope  of  the  beam  at  px  is  tangent  to  12  and  the 
slope  of  the  beam  at  p2  is  tangent  to  23 .  The  beam  element 
consists  of  two  truss  sub-elements  (Pl2  and  2 p2)  and  a 
torsional-spring  bending  sub-element  ( pxlp2 )•  The  internal 
force  in  a  sub-truss  element  is  given  by  equation  (9).  The 
internal  moment  in  the  bending  sub-element  is  given  by: 

M  =  —  Aa  +  —cc  (10) 

4  4 

where  /  is  the  cross-sectional  effective  moment  of  inertia,  L0 
is  the  total  un-stretched  length  of  the  bending  element  which 
is  equal  to  the  length  of  Pl2  plus  2  p2 ,  and  A  a  is  the  change 
in  angle  between  Pl2  and  2Pl  from  the  unstressed 
configuration.  Figure  2b  shows  how  a  beam  is  discretized 
using  the  3-noded  beam  element.  This  thin  beam  element 
does  not  have  a  torsional  response  along  the  axis  of  the 
beam.  In  addition,  it  assumes  that  the  bending  moments  of 
inertia  of  the  cross-section  around  two  perpendicular  cross- 
section  axes  are  the  same. 


Beam  elements 

Figure  2.  (a)  3-noded  beam  element;  (b)  finite  element  discretization 
of  a  beam  using  the  3-noded  beam  element. 
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3.3  Lumped  Parameters  Thick  Spatial  Beam 
Element 


A  lumped  parameter  spatial  beam  element  is  used  to  model 
thick  beams.  The  element  has  two  rigid  body  type  nodes. 
Therefore,  each  node  has  3  translational  DOFs  and  the 
rotation  matrix  is  used  to  measure  the  spatial  rotation  of  the 
node.  The  beam  element  local  X-axis  is  in  the  direction  of 
the  line  connecting  the  two  nodes.  The  Y  and  Z  axes  are 
normal  to  the  X-axis  and  to  each  other.  They  are  defined  as 
an  average  of  the  local  Y  and  Z  axes  of  the  two  nodes.  The 
element  models  each  response  component  independently 
using  simpler  elements  (figure  3): 


•  Axial  response  is  modeled  using  a  zero-length  linear 
spring-damper  truss  element  at  the  center  point  between 
the  two  nodes.  The  internal  force  along  the  X-axis  for  the 
truss  element  is  given  by: 


„  EA  ca- 

F  = - dv  + - dv 


(11) 


where  E  is  the  Young’s  modulus,  C  is  the  damping 
modulus,  A  is  the  cross-sectional  area,  /0  is  the  unstretched 
length  of  the  beam  and  dx  is  the  deflection  between  two 
nodes  along  the  beam’s  X-axis. 


•  Shear  response  is  modeled  using  two  zero-length  linear 
spring-damper  elements  located  at  the  center  point 
between  the  two  nodes.  One  element  produces  a  reaction 
force  along  the  element  7-axis  and  the  other  element 
produces  a  reaction  force  along  the  element  Z-axis.  The 
internal  forces  for  the  7-truss  element  along  the  7-axis  and 
for  the  Z-truss  element  along  the  Z-axis  are  given  by: 


F  =  K 


fGA  ,  CA  ■,  ' 

- u  H - d 

1  y  i  y 
V  *o  *o 


(12a) 


F '  =  K. 


r  GA  ,  CA  ■  ' 

- d  H - d 

1  1  ‘ 

V  *o 


(12b) 


where  G  is  the  shear  modulus,  Ky  and  Kz  are  the  cross- 
section  shear  factors,  and  dy  and  dz  are  the  deflections 
between  two  nodes  along  the  beam’s  7-axis  and  Z-axis, 
respectively. 

•  Bending  response  is  modeled  using  two  torsional  spring- 
damper  elements  located  at  the  center  point  between  the 
two  nodes.  One  element  produces  the  bending  response 
around  the  7-axis  and  the  other  element  produces  the 
bending  response  around  the  Z-axis.  The  internal  moments 
for  the  7  torsional  element  and  for  the  Z  torsional  element 
are  given  by: 


EI  Cl  . 
My=—ay+  —  ay 


*0 


*0 


M  _  =F^a. 


(13a) 

(13b) 


where  Iyy  and  Izz  are  the  second  moments  of  inertia  of  the 
cross-section  around  the  7  and  Z  axes,  respectively,  0Cy  and 
az  are  the  change  in  angles  between  two  nodes  around  the 
beam’s  7-axis  and  Z-axis,  respectively. 

•  Torsional  response  is  modeled  using  a  torsional  spring- 
damper  element  located  at  the  center  point  between  the 
two  nodes.  The  element  produces  the  torsional  response 
around  the  X-axis.  The  internal  moments  for  the  X 
torsional  element  is  given  by: 

,  j  EJ  CJ  .  /i  a\ 


where  J  is  the  torsional  moment  of  inertia  of  the  cross- 
section,  ax  is  the  change  in  angle  between  two  nodes 
along  the  beam’s  X— axis. 

The  deflections  dx,  dy  and  dz  can  be  obtained  as  follows. 
The  center  point  the  element  for  each  of  the  element  nodes 
found  ( Xc\t ,  Xc2 . ,  i  =  1,2,3).  Then,  the  vector  connecting 
the  two  points  is  found: 

Ui2i=Xc2i-Xcii  (15) 


The  dot  product  vector  U 12  with  the  X-axis,  7-axis  and  Z- 
axis  vectors  of  the  beam  give  dx ,  dy  and  dz ,  respectively.  The 
angle  changes  ax,  0Cy  and  az  can  be  found  using  the  nodes 
rotation  matrices. 

The  main  features  of  the  lumped  parameters  beam  element 
are: 

•  It  can  model  spatial  beams  including  bending  along  the 
two  transverse  cross-section  directions  and  torsion. 

•  Shear  deformation  is  modeled. 

•  Since  the  beam  is  modeled  using  nodes  having  both 
translational  and  rotational  DOFs  (therefore,  they  can  be 
thought  of  as  rigid  bodies),  therefore,  detailed  3D  contact 
surfaces  for  the  beam  can  be  used.  This  allows  the  use  of 
this  beam  element  to  model  the  detailed  geometry  of  track 
segments. 
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Figure  3.  Spatial  lumped  parameters  2-node  beam  element.  Each 
node  is  modeled  as  a  rigid  body.  The  X,  Y,  Z  local  axes  of  each  node 
are  shown  red,  green  and  blue,  respectively.  The  element  consists  of 
3  truss  elements  (for  modeling  the  axial  and  shear  response)  and  3 
torsional  spring  elements  (for  modeling  the  bending  and  torsion 
response)  located  at  the  center  point  of  the  element. 

3.4  Natural-Modes  Brick  Element 

The  8-nodes  natural-modes  brick  element  developed  in  [9, 
10]  is  strategically  designed  to  model  all  the  deformation 
and  rigid  body  modes  (figure  4)  of  a  brick  element  while 
avoiding  locking  and  spurious  modes.  All  the  element  nodes 
are  point  mass  type  nodes  with  only  translational  DOFs. 

Two  types  of  sub-elements  are  used  to  model  the  brick: 
two-node  truss  element  and  four-node  surface  shear  element. 
Twelve  truss  elements  along  the  twelve  edges  of  the  element 
are  used  to  model  the  membrane  and  bending  modes  of  the 
element.  Six  surface  shear  elements  are  introduced  at  each  of 
the  six  element  surfaces  to  model  the  shear  and  warping 
modes  of  the  element  (figure  5).  The  derivations  of  the 
stiffness  characteristics  and  structural  forces  generated  by 
the  truss  and  the  surface  shear  elements  are  given  in  .  The 
brick  element  with  the  following  characteristics: 

•  The  deformation  element  modes,  including,  membrane, 
bending,  shear,  and  warping  are  accurately  modeled.  Thus 
the  element  does  not  have  any  spurious  modes. 

•  The  shear  stresses  are  evaluated  as  the  average  stresses 
over  each  of  the  element  faces,  thus  the  element  does  not 
exhibit  shear  locking. 


•  All  element  structural  forces  are  aligned  with  the  element 
edges,  thus  the  element  does  not  exhibit  membrane  or 
membrane-warping  locking. 


3  rotation  rigid  body 
modes 


3  translation  rigid  body 
modes 


3  membrane  (stretch) 
modes 


6  bending  modes 


3  asymmetric  bending 
modes 


3  warping  modes 

Figure  4.  24  rigid  body  and  deformation  modes  of  a  spatial  8-node 
brick  element. 


Twelve  2-node  truss  elements 

Membrane  and  bending  modes  are 
modeled  using  tmss  sub-elements 
along  each  of  the  twelve  element 
edges. 


Six  4-node  surface  shear  elements 

Shear  and  warping  modes  are 
evaluated  using  surface  shear  sub¬ 
elements  on  each  of  the  six  element 
surfaces. 


Figure  5.  Sub-elements  of  the  8-node  lumped-parameters  brick 
element. 


4.  CONTACT  MODEL 

The  penalty  technique  is  used  to  impose  the  normal  contact 
constraints  between  finite  element  nodes  or  points  on  a  rigid 
body  and  finite  element  surfaces  or  quadrilateral  surfaces  of 
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rigid  bodies  [5,  6].  The  first  step  is  to  find  the  position  and 
velocity  of  the  contact  nodes  and  points.  For  finite  element 
nodes  the  global  position  xgp  and  velocity  xgp  of  a  contact 
node  relative  to  the  global  inertial  frame  are  readily 
available: 


XGpi  ~  XKi 

(16a) 

XGP[  =  XKi 

(16b) 

where  xKi  and  xKi  are  the  position  and  velocity  vectors  of 
contact  node  K.  For  rigid  bodies  the  global  position  xgp  and 
velocity  xgp  of  a  contact  point  are  given  by: 

xgp.  —  Xbf.  +  Rbf..  xlp  . 

1  1  LJ  J 

(17a) 

XGPi  —  Xbf1  +  Rbf^  (WbfXXLj?)  . 

(17b) 

where  Xbf  and  Xbf  are  the  global  position  and  velocity 
vectors  of  the  rigid  body’s  frame,  Rbf  is  the  rotation  matrix 
of  the  rigid  body  relative  to  the  global  reference  frame,  Wbf 
is  the  rigid  body’s  angular  velocity  vector  relative  to  its  local 
frame,  and  xlp  is  the  position  of  the  contact  point  relative  to 
the  rigid  body’s  frame. 

The  frictional  contact  force  Fc  at  each  contact  point/node 
(sum  of  the  normal  contact  and  tangential  friction  forces)  is 
transferred  as  a  force  and  a  moment  to  the  center  of  the  rigid 
body.  The  negative  of  this  force  is  transferred  to  the  contact 
surface  element  by  distributing  it  to  the  nodes  forming  the 
surface  using  the  element  shape  function: 

Fu=-NkFCi  (18) 

where  Nk  are  the  surface  element  shape  functions  at  the 

contact  point  and  Fki  are  the  contact  forces  on  node  k  of  the 
surface  element.  In  case  the  contact  body  is  a  rigid  body, 
then  this  force  can  also  be  transferred  to  the  center  of  the 
contacting  rigid  body  as  a  force  and  moment: 

Fi=-FCi  (19a) 

Mi=-(xLPixRBFjiFci )  (19b) 

xlPj  —  Rbf (xGPi  -  Xbf1  )  (20) 

where  Ft  is  the  contact  force  at  the  CG  of  the  contact  rigid 
body  (center  of  the  body  frame),  Mt  in  the  contact  moment 
on  the  contact  rigid  body,  xlcp  is  the  position  of  the  contact 
point  relative  to  the  rigid  body’s  frame  and  xgcP  is  the 
position  of  the  contact  point  relative  to  the  global  reference 
frame.  Thus,  the  contact  algorithm  supports  contact  flexible- 
flexible,  rigid-rigid  and  rigid-flexible  body  contact. 


4.1  Penalty  Normal  Contact  Model 

The  penalty  technique  is  used  for  imposing  the  constraints 
in  which  a  normal  reaction  force  ( Fnormai )  is  generated  when  a 
node  penetrates  in  a  contact  body  whose  magnitude  is 
proportional  to  the  penetration  distance.  In  the  present 
formulation,  the  force  is  given  by  [5,  6]: 

\  Cpd  d>  0 

Fnormai  —  AkpCl  +  A\  .  (21) 

[SpCpd  d  <  0 

d  =  Vn-  ni  (22) 


Contact 


where  A  is  the  area  of  the  rectangle  associated  with  the 
contact  point,  kP  and  cP  are  the  penalty  stiffness  and  damping 
coefficient  per  unit  area;  d  is  the  closest  distance  between 

the  node  and  the  contact  surface  (figure  6);  d  is  the  signed 
time  rate  of  change  of  d;  sP  is  a  separation  damping  factor 
between  0  and  1  which  determines  the  amount  of  sticking 
between  the  contact  node  and  the  contact  surface  at  the  node 
(leaving  the  body);  n  is  the  normal  to  the  surface  and  v«.  is 

the  velocity  vector  in  the  direction  of  n  .  The  normal  contact 
force  vector  is  given  by: 

Fn i  —  Yl ■  Fnormai  (23) 

The  total  force  on  the  node  generated  due  to  the  frictional 
contact  between  the  point  and  surface  is  given  by: 

Fpo  int =  Ft i  +  Fn-  (24) 

4.2  Asperity  Friction  Model 

An  asperity-spring  friction  model  is  used  to  model  joint 
and  contact  friction  [20]  in  which  friction  is  modeled  using  a 
piece-wise  linear  velocity-dependent  approximate  Coulomb 
friction  element  in  parallel  with  a  variable  anchor  point 
spring.  The  model  approximates  asperity  friction  where 
friction  forces  between  two  rough  surfaces  in  contact  arise 
due  to  the  interaction  of  the  surface  asperities.  Ftt  is  the 
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tangential  friction  contact  force  vector  transmitted  to  the 
contact  body  at  the  contact  point.  It  is  given  by: 


search.  At  the  initialization  of  the  algorithm  the  following 
steps  are  performed: 


F tf  —  F tangent  tf 


(25) 


Simple  approximate 
Coulomb  friction  element 


£> 


Spring  with  a 
variable  anchor 
point. 

Figure  7.  Asperity  spring  friction  model. 


© 


F tangent 


The  asperity  friction  model  is  used  along  with  the  normal 
force  to  calculate  the  tangential  friction  force  {F tangent)  [7]. 
When  two  surfaces  are  in  static  (stick)  contact,  the  surface 
asperities  act  like  tangential  springs.  When  a  tangential  force 
is  applied,  the  springs  elastically  deform  and  pull  the 
surfaces  to  their  original  position.  If  the  tangential  force  is 
large  enough,  the  surface  asperities  yield  (i.e.  the  springs 
break)  allowing  sliding  to  occur  between  the  two  surfaces. 
The  breakaway  force  is  proportional  to  the  normal  contact 
pressure.  In  addition,  when  the  two  surfaces  are  sliding  past 
each  other,  the  asperities  provide  resistance  to  the  motion 
that  is  a  function  of  the  sliding  velocity  and  acceleration,  and 
the  normal  contact  pressure.  Figure  7  shows  a  schematic 
diagram  of  the  asperity  friction  model.  It  is  composed  of  a 
simple  piece-wise  linear  velocity-dependent  approximate 
Coulomb  friction  element  (that  only  includes  two  linear 
segments)  in  parallel  with  a  variable  anchor  point  spring. 

4.3  Contact  Search 

Contact  detection  is  performed  between  contact  points  on  a 
“master  contact  surface”  and  a  polygonal  surface  called  the 
“slave  contact  surface”.  The  contact  points  of  the  master 
contact  surface  can  either  be  point  mass  nodes  or  points  on  a 
contact  surface  of  a  rigid  body  type  node.  The  slave  contact 
surface  can  be  a  polygonal  surface  connecting  point  mass 
type  nodes  or  a  polygonal  surface  on  a  rigid  body  type  node. 
Contact  between  the  contact  points  of  the  master  surface  and 
the  polygons  of  the  slave  surface  is  detected  using  a  binary 
tree  contact  search  algorithm  which  allows  fast  contact 


•  Each  slave  polygonal  contact  surface  is  divided  into  2 
blocks  of  polygons.  The  bounding  box  for  each  block  of 
polygons  is  found.  Then  each  of  those  blocks  of  polygons 
is  divided  into  2  blocks  and  again  the  bounding  boxes  for 
those  blocks  are  found.  This  recursive  division  continues 
until  there  is  only  one  polygon  in  a  box. 

•  For  each  master  contact  surface  the  contact  points  are 
divided  into  2  blocks.  The  bounding  sphere  for  each  block 
of  points  is  found.  Then,  each  of  those  blocks  of  points  is 
divided  into  2  blocks  and  again  the  bounding  spheres  for 
those  blocks  are  found.  This  recursive  division  continues 
until  there  is  only  one  point  (with  a  bounding  sphere  of 
radius  0). 

During  the  solution  the  following  steps  are  performed.  For 
each  master  contact  sphere,  the  radius  of  the  contact  sphere 
is  added  to  the  size  of  the  bounding  box,  and  then  we  check 
if  the  center  point  of  the  sphere  is  inside  a  bounding  box.  If 
the  center  of  the  contact  sphere  is  not  inside  any  bounding 
box,  then  all  the  points  inside  that  sphere  are  not  in  contact 
with  the  surface.  If  the  center  of  the  contact  sphere  is  inside 
a  bounding  box  then  the  two  sub-bounding  boxes  are 
checked  to  determine  if  the  point  is  inside  either  one.  If  it  is, 
then  the  sub-contact  spheres  are  checked.  If  a  contact  point 
is  found  to  be  inside  the  lowest  level  bounding  box,  then  a 
more  computationally  intensive  contact  algorithm  between  a 
point  and  a  polygon  is  used  to  determine  the  depth  of  contact 
and  the  local  position  of  the  contact  point  on  the  polygon. 

This  search  algorithm  has  a  theoretical  average 
computational  cost  of  log(m)  x  log (n),  where  m  is  the 
number  of  points  of  the  master  surface  and  n  is  the  number 
of  polygons  of  the  slave  surface.  It  allows  detecting  contact 
between  surfaces  containing  millions  of  polygons  and 
contact  points  in  near  real-time. 

5.  JOINT  CONSTRAINTS 

Each  rigid  body  can  have  a  number  of  connection  points. 
A  connection  point  is  a  point  on  the  body  where  joints  can 
be  located.  A  connection  point  does  not  add  additional  DOFs 
to  the  system  .The  connection  point  can  be: 

•  A  point  mass  type  node. 

•  A  point  on  a  rigid  body. 

•  An  arbitrary  point  inside  a  finite  element. 

5.1  Connection  point  location 

If  the  connection  point  is  a  node  then  equations  (16a)  and 
(16b)  are  used  to  find  the  global  position  xgp  and  velocity 
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x gp  of  the  connection  point.  If  the  connection  point  is  a 
fixed  point  on  rigid  body  B  then  equations  (17a)  and  (17b) 
are  used  to  find  the  global  position  and  velocity  of  the 
connection  point.  If  the  connection  point  is  a  point  inside  a 
finite  element,  then  xgp  and  velocity  xgp  are  given  by: 

(26) 

xcp'^NAMx'j,  (27) 

where  J  is  the  local  node  number  of  the  element,  Nj  (^ )  are 
the  interpolation  functions  of  the  element,  %  are  the  natural 
element  coordinates  of  the  fixed  point  and  Xj.  is  the  position 

vector  of  local  node  J  of  the  element  relative  to  the  global 
reference  frame. 

5.2  Spherical  joint  constraint  force 

A  joint  is  defined  by  defining  the  relation  between 
connection  points.  For  example,  a  spherical  joint  between 
two  connection  points  is  defined  as: 

xc^xd]  (28) 

where  xc\\  is  the  position  vector  of  the  first  point  and  xd]  is 


the  position  vector  of  the  second  point, 
imposed  using  the  penalty  technique  as: 

This  constraint  is 

Fc  =  kP\v\+cP  v.v. 

(29) 

V-  =  Xc\\  -  Xc2 \ 

(30) 

V-  =  Xcl]  ~  Xcl\ 

(31) 

FC{  =  Fc  v. 

(32) 

where  FCi  is  the  penalty  reaction  force  on  the  connection 
point,  kP  is  the  penalty  spring  stiffness,  and  cP  is  the  penalty 
damping.  The  constraint  force  is  applied  on  the  two 
connection  points  in  opposite  directions.  Depending  on  the 
type  of  connection  point  the  constraint  force  is  applied  as 
follows.  If  the  connection  point  is  a  point  on  a  rigid  body, 
then  it  is  transferred  to  the  node  at  the  center  of  the  body  as  a 
force  and  a  moment  using  equations  (14a)  and  (14b).  If  the 
connection  point  is  a  node,  then  the  constraint  force  is 
applied  directly  to  the  node: 

FL=Fc‘  (33) 

If  the  connection  point  is  a  point  inside  a  finite  element,  then 
it  is  applied  to  the  nodes  of  the  element  using: 

Fji  =  Nj{^j)  Fc\  (34) 


where  J  is  the  local  node  number  of  the  element,  is  the 

force  on  local  node  J  of  the  element  relative  to  the  global 
reference  frame. 

Using  equations  (16,  17,  and  26-34),  the  following  types 
pin  joints  can  be  modeled: 

•  Spherical  -joint  between  two  rigid  bodies. 

•  Spherical-joint  between  a  rigid  body  and  a  finite  element 
point. 

•  Spherical-joint  between  a  rigid  body  and  a  point  particle 
type  node. 

•  Spherical-joint  between  two  element  points. 

•  Spherical-joint  between  an  element  point  and  a  point 
particle  type  node. 

•  Spherical-joint  between  two  nodes. 

The  constraint  forces  are  applied  to  the  connection  point 
node(s)  by  assembling  them  into  the  global  structural  forces 
Fs  in  equation  (1).  Also,  the  constraint  moments  are  applied 
to  the  nodes  by  assembling  them  into  the  global  structural 
torques  Ts  in  equation  (2). 

Revolute  joints  can  be  modeled  by  placing  two  spherical 
joints  along  a  line.  Other  types  of  joints  such  as  prismatic, 
cylindrical,  universal  and  planar  joints  can  also  be  modeled 
by  writing  the  constraint  equation,  then  writing  the 
corresponding  penalty  forces  and  moments  on  the 
connection  points. 

6.  EXPLICIT  SOLUTION  PROCEDURE 

The  solution  fields  for  modeling  multibody  systems  are 
defined  at  the  model  nodes.  Note  that  a  rigid  body  modeled 
as  one  finite  element  node.  These  solutions  fields  are: 

•  Translational  positions. 

•  Translational  velocities. 

•  Translational  accelerations. 

•  Rotation  matrices. 

•  Rotational  velocities. 

•  Rotational  accelerations. 

The  explicit  time  integration  solution  procedure  predicts 
the  time  evolution  of  the  above  response  quantities.  An 
advantage  of  explicit  solution  procedures  is  that  they  are 
“embarrassingly”  parallel.  The  procedure  described  below 
achieves  near  linear  speed-up  with  the  number  of  processors 
on  shared  memory  parallel  computers.  The  procedure  is 
implemented  in  a  multibody  dynamics  finite  element  code 
and  is  outlined  below: 

1  -  Prepare  the  run: 

a.  Set  the  initial  conditions  for  the  solution  fields 
identified  above. 
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b.  Create  a  list  of  all  the  finite  elements  (Those  also 
include  joints  and  master  contact  surfaces  which  are 
considered  elements). 

c.  Create  a  list  of  elements  that  will  run  on  each 
processor.  This  is  done  using  an  algorithm  which  tries 
to  make  the  computational  cost  on  each  processor 
equal. 

d.  Create  a  list  of  all  the  prescribed  motion  constraints. 

e.  Calculate  the  solid  masses  for  each  finite  element  node 
by  looping  through  the  list  of  finite  elements.  Note  that 
the  masses  are  fixed  in  time. 

f.  Loop  over  all  the  elements  and  find  the  minimum  time 
step  for  the  explicit  solution  procedure. 

2-  Loop  over  the  solution  time  and  increment  the  time  by 
At  each  step  while  doing  the  following: 

a.  Set  the  nodal  values  at  the  last  time  step  to  be  equal  to 
the  current  nodal  values  for  all  solution  fields. 

b.  Do  2  iterations  (a  predictor  iteration  and  a  corrector 
iteration)  of  the  following: 

i.  Initialize  the  nodal  forces  and  moments  to  zero. 

ii.  Calculate  the  nodal  forces  and  moments  by  looping 
through  all  the  elements  (and  joints)  while 
calculating  and  assembling  the  element  nodal 
forces.  This  is  the  most  computational  intensive 
step.  This  step  is  done  in  parallel  by  running  each 
list  of  elements  identified  in  step  l.c  on  one 
processor. 

iii.  Find  the  nodal  values  at  the  current  time  step  using 
the  semi-discrete  equations  of  motion  and  the 
trapezoidal  time  integration  rule  (Eqs.  1-5). 

iv.  Execute  the  prescribed  motion  constraints  which  set 
the  nodal  value(s)  to  prescribed  values. 

v.  Go  to  the  beginning  of  step  2. 

7.  EXPERIMENTAL  SETUP 

The  model  is  validated  using  a  full-scale  army  tracked 
vehicle.  The  track  was  removed  and  the  vehicle  was  placed 
on  an  n-post  motion  base  simulator  in  the  US  army’s  Tank 
Automotive  Research,  Development  and  Engineering 
Center4 s  (TARDEC)  Simulation  Laboratory  (TSL)  (see 
figures  1  and  8).  The  n-post  motion  simulator  consists  of 
linear  hydraulic  actuators  each  placed  under  one  of  the 
vehicle  road  wheels.  The  road- wheels  are  connected  to  the 
road  arms  using  re  volute  joints.  The  vehicle  suspension 
system  consists  of  the  road  arms  and  torsion  bars  between 
the  roads  arms  and  chassis.  The  road  arms  are  connected  to 
the  chassis  using  revolute  joints.  Each  road  arm  is  free  to 
rotate  between  two  hard  suspension  stops  at  0  degrees  (when 
the  arm  is  horizontal)  and  about  60  degrees  from  the 
horizontal.  Each  actuator  has  one  degree  of  freedom  along 
the  vertical  direction  and  can  be  independently  commanded 


to  follow  a  certain  vertical  displacement  time-history.  Thus, 
the  actuators  can  be  controlled  in  such  a  way  as  to  simulate 
vertical  position  time-histories  of  the  wheels  as  the  vehicle 
goes  over  a  rough  road  that  has  bumps  and  potholes.  Note 
that  the  motion  simulator  cannot  simulate  the  inertial 
centrifugal  inertia  forces  that  arise  due  to  the  time-varying 
motion  of  the  vehicle  on  the  road. 


Figure  8.  Snapshot  of  the  computational  model  of  the  tracked 
vehicle  on  the  n-post  motion  simulator. 


Figure  9.  Dishpans  with  front  and  back  cylindrical  stops  and  side 
guard  rails. 

The  vehicle  in  our  experimental  setup  has  12  road  wheels, 
6  on  each  side.  We  number  the  wheels  from  front  to  back  1 
through  6.  Each  hydraulic  actuator  is  composed  of  a 
cylinder,  piston  and  a  dishpan  on  which  the  wheel  rests.  The 
dishpans  for  wheels  1,  3,  4  and  6  (colored  silver)  have  front 
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and  back  horizontal  cylindrical  stops  as  well  as  side  guard 
rails  (figure  9).  Those  dishpans  are  used  in  order  to  ensure 
that  the  vehicle  does  not  “slide  off’  the  motion  simulator. 
Flat  dishpans  are  used  for  wheels  2  and  5  (figure  10).  In 
addition,  a  strap  harness  was  attached  to  the  top  of  the 
vehicle  during  the  experiments  for  safety  reasons.  The 
harness  was  loose  so  as  not  to  interfere  with  the  motion  of 
the  vehicle. 


Figure  10.  Flat  dishpans. 


The  angles  of  the  12  roads  arms  are  measured  and  sampled 
at  a  rate  of  205  samples/sec.  The  experiment  12  road  arms 
angles  measurements  are  compared  with  the  results 
predicted  using  the  computational  model.  A  camera  mounted 
on  the  ground  is  used  to  record  the  motion  of  the  vehicle. 
The  camera  is  set  to  capture  30  frames/sec.  The  camera  view 
is  shown  in  figure  8. 

The  torsion  bars  at  the  road  wheels  were  charged  such  that 
the  road  arms  are  initially  half  way  between  the  two 
suspension  stops  (i.e.  at  about  30  degrees  from  the 
horizontal)  when  the  vehicle  is  at  rest  under  gravity  and  the 
track  is  off.  Figure  1 1  shows  the  rotational  stiffness  (torque 
versus  angle)  of  the  torsion  bar  at  each  road  wheel.  The 
torsion  bars  rotational  stiffness  was  experimentally  obtained. 

The  matrix  of  excitations  applied  to  the  vehicle  is  shown  in 
Table  1.  Two  main  types  of  excitations  are  used:  shake  and 
pitch.  In  the  shake  excitation  all  the  actuators  are  moving  the 
same  way  (i.e.  their  motions  are  in  phase).  In  the  pitch 


excitation  the  front  and  rear  actuators  are  180°  out  of  phase. 
The  intermediate  actuators  are  moved  in  such  as  way  that  at, 
at  point  in  time,  the  line  connecting  the  center  of  the  front 
dishpan  and  center  of  the  rear  dishpan  passes  through  the 
center  of  the  intermediate  dishpans.  This  way  the  vehicle  is 
rotating  around  its  pitch  axis.  The  pitch  and  shake 
excitations  simulate  the  vehicle  going  over  bumps  and 
potholes.  For  each  motion  type,  harmonic  excitations  at 
amplitudes  40,  60  or  70  mm  and  frequencies  from  0.6  to  1.5 
Hz  are  applied  to  the  vehicle.  Each  excitation  is  applied  for 
about  1 5  sec  with  the  initial  3  sec  and  the  last  3  sec  used  to 
ramp  up  and  down  the  excitation  amplitude.  Also,  for  each 
ramp  motion  type,  a  ramp  excitation  of  90  mm  is  applied  to 
the  vehicle,  where  the  vehicle  is  raised  90  mm  in  0.2  sec, 
then  the  actuators  dwell  for  5  sec,  finally  the  vehicle  is 
lowered  90  mm  in  0.2  sec. 


Figure  11.  Experimetally  measured  rotational  stiffness  (torque 
versus  angle)  of  the  torsion  bar  at  each  road  wheel. 


Table  1.  Matrix  of  the  tracked  vehicle  experiments  (38  experiments). 


Excitation 

Amplitude 

Frequency  (Hz) 

Ramp 

Duration  (sec) 

Harmonic  - 
Pitch 

40  mm 

0.6;  0.7;  0.8;  0.9;1.0;  1.1;  1.2;  1.3;  1.4 

- 

60  mm 

1.1;  1.2 

- 

70  mm 

0.6;  0.7;  0.8;  0.9;1.0 

- 

Harmonic- 

Shake 

40  mm 

0.6;  0.7;  0.8;  0.9;  1.0;  1.1;  1.2;  1.3;  1.4;  1.5 

- 

70  mm 

0.6;  0.7;  0.8;  0.9;  1.0;  1.1;  1.2;  1.3;  1.4;  1.5 

- 

Ramp  - 
Pitch 

90  mm 

- 

0.2 

Ramp  - 
Shake 

90  mm 

- 

0.2 
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8.  VALIDATION  STUDY 

The  experiments  in  Table  1  were  simulated  using  the 
present  model.  First  the  experiments  were  used  to  find  an 
average  damping  moment  versus  angular  velocity  for  the 
suspension  system’s  torsion  bars  (figure  12).  This  average 
damping  moment  versus  angular  velocity  curve  was  used  for 
all  the  torsion  bars.  The  damping  moment  versus  angular 
velocity  typically  depends  on  the  oil  charge  of  the  torsion 
bar.  Note  that  the  curve  is  nonlinear.  For  very  small  angular 
velocities,  the  tension  bar  produces  a  large  amount  of 
damping.  For  larger  values  of  angular  velocity,  the  damping 
moment  is  nearly  linear  with  angular  velocity,  and  gradually 
deviates  from  linearity  and  becomes  stiffer  for  very  large 
angular  velocities.  Also,  from  the  experiments  we  found  that 
there  is  small  friction  moment  in  the  tension  bars  of  about  5 
N.m.  Note  that  the  present  model  accounts  for  non-linear 
stiffness,  damping  and  friction  of  the  suspension  system. 
Also,  the  model  allows  the  wheels  to  leave  the  dishpans. 


Figure  12.  Damping  moment  versus  angular  velocity  for  the  torsion 

bars. 


The  graphs  in  Figures  13-20  show  the  results  of  4 
representative  runs  from  the  set  of  38  runs  in  the  test  matrix 
shown  in  Table  1.  Those  figures  show  comparisons  of  the 
road  arm  angles  between  the  experiment  and  the  DIS 
simulation.  The  average  difference  in  magnitude  between 
the  experiment  and  simulation  is  about  15%.  The  average 
response  difference  between  the  right  and  left  sides  of  the 
vehicle  was  also  about  15%.  The  main  sources  of  error 
between  the  experiments  and  simulation  in  order  of 
importance  are: 

•  At  the  higher  frequency  shake  excitations  (1.4,  1.5  Hz) 
and  pitch  excitations  (1.0,  1.1  and  1.2  Hz),  the 


amplitude  of  motion  of  the  road  arms  exceeded  ±6°.  At 
those  amplitudes,  the  wheels  where  bumping  and 
“climbing  over”  the  dishpans  cylindrical  stops.  The 
experimental  results  for  those  runs  cannot  be  accurately 
modeled  because  the  position  of  a  wheel  relative  to  the 
dishpans  stops  is  not  known  and  varies  between  runs. 
For  those  excitations  the  difference  between  the 
experiment  and  simulation  was  on  the  average  about 
30%  with  a  maximum  difference  of  about  50%.  Also, 
note  that  the  range  of  motion  of  a  road  arm  is  about 
±29°,  but  the  highest  motion  amplitude  that  we  tested 
was  about  ±10  degrees.  We  could  not  excite  the  vehicle 
beyond  ±10  degrees  road  arm  deflections  for  safety 
reasons  since  the  wheels  started  bumping  and  going 
over  the  stops  at  ±6°. 

•  Asymmetry  of  the  stiffness  and  damping  between  the 
two  sides  of  the  vehicle.  The  difference  between  the 
response  at  the  two  sides  of  the  vehicle  was  on  average 
about  8%. 

•  Presence  of  the  strap  harness. 


9.  TYPICAL  SIMULATIONS 

The  present  model  can  be  used  to  predict  the  dynamic 
response  of  tracked  vehicle.  This  includes  predicting  the 
following: 

•  Stability  of  the  vehicle  while  moving  at  a  certain  speed 
over  a  rough  terrain  or  while  steering. 

•  Maximum  positive  or  negative  obstacles  that  the  vehicle 
can  overcome.  This  includes  initial  speed  and  engine 
torque  required  to  overcome  the  obstacle. 

•  Internal  forces  in  the  various  components  of  the  vehicle 
including  the  track  while  going  over  various  terrains  and 
obstacles.  Those  forces  can  then  be  applied  to  finer 
finite  element  models  of  the  components  for  detailed 
stress  analysis.  The  predicted  stresses  can  be  used  to 
determine  the  whether  the  component  will  fail  under  the 
loads  as  well  as  the  fatigue  life  of  the  component. 

•  Vehicle  energy  and  power  requirements  for  going  at  a 
prescribed  speed  over  a  known  course. 

Further  experimental  tests  are  needed  to  validate  and  tune 
the  model  parameters  to  accurately  predict  the  above 
response  quantities.  Those  experiments  should  include  actual 
road  tests  using  an  instrumented  tracked  vehicle. 
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Figure  14.  Experiment  and  simulation  angles  of  the  road  arms  versus  time  for  input  motion:  harmonic  shake  1.0  Hz  70  mm. 
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Figure  16.  Experiment  and  simulation  angles  of  the  road  arms  versus  time  for  input  motion:  harmonic  shake  1.5  Hz  70  mm. 
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Figure  18.  Experiment  and  simulation  angles  of  the  road  arms  versus  time  for  input  motion:  harmonic  pitch  0.8  Hz  70  mm. 
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Figure  19.  Input  motion:  ramp  shake  90  mm. 


Figure  20.  Experiment  and  simulation  angles  of  the  road  arms  versus  time  for  input  motion:  ramp  shake  90  mm. 
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Figure  21.  Snapshots  of  the  motion  of  the  tracked  vehicle  with  a 
continuous  belt  track.  The  track  is  modeled  using  brick  elements  for 
modeling  the  belt  rubber  and  thin  beam  elements  for  modeling  the 
reinforcements,  (a)  vehicle  going  over.semi-circular  bumps;  (b) 
vehicle  going  over  a  rough  course. 


15 


Figure  22.  Snapshot  of  the  motion  of  the  tracked  vehicle  with  a 
continuous  belt  track  modeled  using  thick  beam  elements. 


Finite  element  models  of  a  typical  tracked  vehicle,  with  a 
continuous  belt-type  track,  are  shown  in  Figures  21  and  22. 
In  those  models  the  chassis,  road  arms,  road  wheels,  idlers 
and  sprockets  are  modeled  using  rigid  bodies.  The 
suspension  rotary  spring-dampers  are  modeled  as  torsional 
springs.  Rotary  actuators  along  with  control  laws  are  used  to 
model  the  engine  and  steering  system  of  the  vehicle.  In 
Figure  21  the  track  is  modeled  using  brick  elements  for  the 
rubber  matrix  and  thin  beam  for  the  reinforcements.  In 
Figure  22  the  track  is  modeled  using  the  thick  beam 
elements  presented  in  Section  3.3.  Figures  21a  shows  a 
snapshot  of  the  tracked  vehicle  going  over  semi-circular 
bumps.  Figure  21b  shows  a  snap  shot  of  the  tracked  vehicle 
going  over  a  rough  course.  Figure  22  shows  a  snapshot  of 
the  tracked  vehicle  undergoing  a  lane-change. 

10.  CONCLUDING  REMARKS 

A  finite  element  model  for  predicting  the  fully  coupled 
dynamic  response  of  flexible  multibody  systems  for  tracked 
vehicles  was  presented.  The  model  has  the  following 
characteristics: 

•  Parallel  explicit  time-integration  solver. 

•  Algorithm  for  accurate  accounting  for  the  rigid  body 
rotational  motion.  The  rigid  bodies  rotational  equations 
of  motion  are  written  in  a  body-fixed  frame  with  the  total 
rigid  body  rotation  matrix  updated  each  time  step  using 
incremental  rotations. 

•  Total-Lagrangian  lumped  parameters  3D  finite  elements 
including  spring/truss,  thin  beam,  thick  beam,  thick  shell, 
and  solid  finite  elements. 

•  Penalty  technique  for  modeling  joint  constraints 
including  spherical,  revolute,  cylindrical  and  prismatic 
joints. 

•  Penalty  technique  for  modeling  normal  contact. 

•  An  asperity  friction  model  is  used  to  model  joint  and 
contact  friction. 

•  General  fast  hierarchical  bounding  box-bounding  sphere 
contact  search  algorithm  for  finding  the  contact 
penetration  between  points  on  master  contact  surfaces 
and  polygons  on  slave  contact  surfaces. 

•  Continuous  belt-type  tracks  are  modeled  using  brick 
elements  representing  the  rubber  matrix  and  beam 
elements  along  the  width  and  circumference  of  the 
belt/segment  for  modeling  the  belt  reinforcements. 

•  Continuous  or  segmented  tracks  are  modeled  thick  beam 
elements. 

A  validation  study  of  the  finite  element  model  was  carried 
out  using  a  tracked  vehicle  (without  the  track)  mounted  on 
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an  n-post  motion  simulator.  The  study  shows  that  the  model 
can  predict  reasonably  well,  within  15%  difference  on 
average,  the  response  of  the  vehicle.  The  maximum 
difference  between  the  experiments  and  simulation  was  50% 
for  the  large  frequency  excitations,  mainly  due  to  the  wheels 
climbing  over  the  dishpan  stops  in  the  experiments. 
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